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ABSTRACT 



The Finite Element Analysis Program (FEAP) was expanded 
to solve linear and nonlinear, two and three dimensional 
heat conduction problems. The usual types of . boundary 
conditions, including radiation, may be specified. A wide 
range of two- and three-time level schemes for the solution 
of time dependent problems is available in the program and a 
discussion of those most commonly used is presented. The 
algorithms for the solution of typical practical problems are 
described and several numerical examples are presented. The 
results are compared with the available analytical 
solutions. A listing of this expanded version of FEAP and 
the corresponding user's instructions are provided. 
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I 



INTRODUCTION 



This thesis describes an application of the Finite 
Element method to Heat Transfer analysis through the use of 
a computer program. This program was designed to solve 
linear and nonlinear, steady and unsteady, two and three 
dimensional heat conduction problems involving temperature 
dependent thermophysical properties and complicated 
radiation/convection boundary conditions. 

The Finite Element Analysis Program (FEAP), programmed 
by Professor R. L. Taylor in the Department of Civil 
Engineering of the University of California, Berkeley, 
served as the point of departure for the present code. It 
has now been expanded with two additional modules for Heat 
Transfer analysis and time integration of first order 
equations, respectively, but the original characteristics of 
a research and educational tool in which the various modules 
can be changed or added to as desired, were maintained. 

Although all equations are retained in core, the program 
can handle realistic engineering problems with several 
hundred unknowns in most computer systems. 

While the algorithms to solve steady heat transfer 
problems are well studied and defined, more research has to 
be done on the solution of unsteady problems. The 
Zienkiewicz two- and three-time level schemes were 
integrated in the program in such a way that a wide range of 
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choices of time Integration algorithms is available to the 
user. A limited study was performed in order to determine 
the characteristics of the most used members of those 
families of numerical schemes. 

The FEAP code is written in FORTRAN IV language and this 
version was constructed and tested on an IBM 360/67 computer 
system with OS/67 release 18. Other systems and 
installations will require modifications. 
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I I . STATEMENT OF OBJECTIVES 



A. EQUATION OF CONDUCTION. INITIAL AND BOUNDARY 
CONDITIONS 



The problem considered for solution 

developed by Arpaci in Ref. 1 and is 

described in the region ft by the equation 

P c |i = V . (EVT) + Q 

subjected to boundary conditions 

T = T on T i 

w 1 

and 



i s thorough 1 y 
mathemat i ca 1 y 

( 1 ) 

( 2 ) 



(£VT)n+q+q c +q r =0 on r 2 (3) 
and initial condition 

lim T = T (4) 
t+0 0 

V is the gradient operator, Ti and r 2 are mutually 
exclusive parts of the boundary of the region ft , T is the 
temperature and t is the time. The thermal capacity pc, the 
thermal conductivity 1< and the rate of internal heat 
generation Q are thermophysical properties dependent on 
temperature. 

In equation (3) rT is an outward unit vector normal to 
the boundary surface while q, q c , q r represent respectively 
the imposed heat flux and the rates of heat flow per unit 
area due to convection and radiation defined as 



q = h 
^c c 



( T - T ) 
ac ' 



and 



( 5 ) 



9 



( 6 ) 



q r ■ FO ( T 4 - T^ r ) - h r C T - T ar ) 

In C 5 ) h c Is the temperature dependent convective heat 

transfer coefficient. In (6) the parameter h is defined by 

r 

the expression 

h r = Fa ( T 2 + l\ x ) ( T + T ar ) , where F is 
the radiant exchange factor and o the Stefan-Bol tzman 
constant. T ac and T ar are the equilibrium temperatures for 
which, respectively, no convection and radiation occurs. 

B. NUMERICAL SOLUTION. DISCRETIZATION BY THE GALERKIN 
METHOD. MATRIX FORMULATION 

The spacewl se discretization of equation (1) using 
cartesian coordinates can be acomplished by Galerkln's 
principle as shown by Zlenkiewicz in Ref. 2 and Lew in Ref. 
3. 

Let the unknown function T be approximated, throughout 

the solution domain at any time t, by the relationship 

T = l N (x,y,z) T. (t) = <N> {T} C 7 ) 

i=l 1 1 

where N^ are the usual shape functions defined piecewise 
element by element, T^ being the nodal parameters. 

The result is 

[K] {T> + [C] {T> + {F} = {0} (8) 

where [K] and [C] are symmetric matrices defined, on the 
element level, as 

W e = J a . 

<6N 
6z 

f f h + h ) <N> <N> dr 
J ^ c r 
l 2 



(<i> <n> - 

> <-^> k ) dfl + 
6z z J 



T 

+ <f> 

6y 



<|N> k 

6y ) 
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( 10 ) 



[C] e = / pc<fJ <N> dfl 






The load vector F is 



{ F } e = -/ <N^ Q dfi + 



ST 

V 

1 2 



( 



h T 
c ac 



h T ) dr 
r ar ' 



( 11 ) 



where the surface integrals are performed only over the 
surfaces where the prescribed boundary condition applies. 

It must be noted that the set of equations (8) is 
nonlinear since the matrices [K] , [C] and vector {F}are 
dependent on T . 

The system of equations (8) may now be solved by any 
automated numerical technique. Code FEAP was written for 
this purpose and Is described in the rest of this thesis. 
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i i i . FINITE ELEMENT ANALYSIS PROGRAM 



The Finite Element Analysis Program (FEAP) was 
programmed and published by Taylor in Ref. 4 where it is 
well explained. The reader should consult this reference for 
detai 1 s about i t. 

This thesis is an expansion of the original program 
towards the solution of Heat Conduction problems and as such 
only the main features of the basic program will be referred 
in this section. A more detailed explanation is reserved for 
the new subroutines added to the program. The complete 
user's instructions for FEAP are given in Appendix A and 
they complement this section. 

A. MAIN PROGRAM 

FEAP can be separated into two basic parts: 

a) Data input module and preprocessor 

b) Solution and output modules 

The data input module must transmit sufficient 
information to the other modules so that each problem can be 
solved. All the input data is stored in a single array 
(integer array M) which is partitioned to store all the data 
arrays# as well as the global arrays# e.g.# stiffness# mass# 
load# etc. 

The total capacity of the program is controlled by the 
dimension of the array M in the blank common of the main 
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program and the corresponding value of the variable MAX. 

The partition of M is performed in the control routine 
(subroutine PCONTR) and the data used for it is supplied by 
the user in the title and control Information cards (first 
two cards for every run of the program). 

Then the nodal coordinates, element connections, 
material properties, nodal loading and boundary conditions 
may be input using the macro control statements. They are 
interpreted by subroutine PMESH and each control is a 
function independent of the others. 

To clarify what was said till this point, the discussion 
of a simple example may be useful. Consider the mesh 
presented in figure 1 and assume that quadratic rectangular 
elements as the one defined in Figure 2 are used. This 
example is completly solved in Appendix B and will be used 
throughout this text. 

The sector of a circular ring was divided in six 
elements and the 33 nodes were numbered. The order of 
numbering is not crucial but in order to improve the profile 
of non-zero coefficients the following general rule should 
be used. 

The numbering should be such as to minimize the nodal 
difference for each element (maximum node nummber minus 
minimum node number). 

The user can now proceed to the preparation of data for 
the program. The first step consists of specifying the 
problem title and the control information. The latter is 

1 Figures are grouped at the end pp 47-79 
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clearly 33 nodes, six elements, one material set, two 
spatial dimensions and eight nodes per element. If this is a 
Heat Transfer problem we have one unknown per node 

.( temperature ) . 

The program expects now the data cards for mesh 
description. Any analysis require at least 

a) coordinate data which follows the macro command COOR 

b) connectivity table which follows the macro command 

ELEM 

c) material data which follows the macro command MATE 

In addition, most analysis will require specification of 
nodal boundary restraint conditions, macro BOUN, and the 
corresponding nodal force or displacement values, macro 

FORC. The term displacement refers to the value of the 
unknown for each specific problem. In a Heat Transfer 

problem the unknown is temperature while in an Elasticity 
analysis it is a geometrical displacement. 

The TEMP macro command (temperature data) shall not be 
used in a Heat Transfer analysis to specify nodal 

temperature. It may be used to specify auxilary nodal 

quant i ties. 

If in our problem we assume the line defined by nodes 
31, 32 and 33 at the temperature of 20, the line defined by 
nodes 1, 2 and 3 at 820 and insulated elsewhere, the macro 
control statements and correspondent data cards will be 
COOR 

1 5.1 
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31 .3 

4 5.116666667 

29 .283333333 



2 


5 .1 


3. 


32 


.3 


3. 


3 


5 .1 


6. 


33 


.3 


6. 


5 


5.116666667 


6. 


30 


.283333333 


6. 


(blank card) 




POLA 






1 


33 1 




(blank card) 




ELEM 






1 


116 


8 


(blank card) 




BOUN 






1 


1 -1 




3 


-1 




31 


1 -1 




33 


-1 




(blank card) 




FORC 






1 


1 320. 




3 


820, 




31 


1 20. 




33 


20. 
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(blank card) 

MATE 

1 2 

10. 10. 8000. 250. 3 

END 

In spite of the fact that the MATE command must be 
included in this part of the data preparation, its 

discussion is left to the section where the element modules 
are treated. 

After completing the mesh data input, we are ready to 
Initiate the problem solution. 

FEAP has modules for variable algorithm capabilities and 
which, if necessary, can be modified or expanded. The basic 
aspect of the variable algorithm program is a macro 
instruction language which can be used to construct modules 
for specific algorithms as needed. This language is 
interpreted by subroutine PMACR and the complete list of 
macro instructions available in this version of FEAP is 
given in Appendix A. 3. 

Here, as an example, we will discuss in detail the 
algorithm for the solution of linear steady state problems 
using the original explanation given in Ref. 4. 

To use the macro programming commands, the user only 
needs to learn the mnemonics of the language. If one wishes 
to form the global stiffness matrix the program instruction 
TANG is used (TANG is the mnemonic for a symmetric tangent 
stiffness matrix and for nonlinear elements would form and 
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assemble into the global stiffness the element tangent 
stiffness computed about the current displacement state; for 
linear elements this is just the linear stiffness matrix). 
For a problem with an unsymmetric tangent stiffness the 
macro command UTAN is used. 

If one wishes to form the right hand side of the 
equations modified for specified displacements one uses the 
program instruction FORM. The resulting equations are 
solved using the Instruction SOLV. 

Printed output can be obtained using the macro command 
DISP, 

The above Instructions are sufficient to solve linear 
steady state problems, that is, the macro instructions 
TANG ( or UTAN ) 

FORM 
SOLV 
Dl SP 

are precisely the required Instructions to solve any linear 
steady state problem. 

If the same block of Instructions is repeated twice or 
more, considerable effort Is wasted in preparing the macro 
Instruction data. To rectify this, looping commands are 
Introduced as the Instruction pair 
LOOP n 

NEXT 

which indicates that looping over all instructions between 
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LOOP and NEXT will occur n times. 

Many other classes of problems can be solved using the 
macro instruction list given in Appendix A. With the above 
short explanation and the discussion presented in section IV 
for Heat Transfer problems, the reader will be able to 
construct his own algorithms using the macro programming 
language. 

B. FINITE ELEMENT SOLUTION MODULES 

Most of FEAP is common to a wide range of problems, but 
some of the computations are linked to the type of problem 
to be sol.ved. It is clear that if one wishes to solve an 
Elasticity problem, the calculation of the stiffness matrix 
will be different from that used for a Heat Transfer 
problem. It Is also clear that differences also arise if one 
Is using triangular elements instead of quadratic 
isoparametric elements as a way of describing the continuum. 

The program was designed such that all computations 
associated with any type of element are contained in an 
element subroutine called ELMTnn where nn is between 01 and 
05 in this version of FEAP. Each element type to be used is 
specified as part of the material property data, following 
macro control MATE. 

Since the objective of this work is the solution of Heat 
Transfer problems, the element modules discussed here are 
the two and three dimensional heat tranfer modules named 
ELMT02 and ELMT03. 
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1. Two Dimensional Heat Transfer Element 



Subroutine ELMT02 is accessed using the number two 
in column ten of each material number card. 

According to the value of the variable I SW, defined 
in subroutine PMACR, a specific function is performed. If 
ISW=1 it reads and prints the material property data and 
line boundary conditions. The stiffnes and mass matrices are 
computed if ISW=3 and I SW=5 respectively. The load vector 
due to internal forces and boundary conditions is calculated 
when ISW=6. 

The material properties may be input as constants or 
as a table of temperature and property values. If any other 
way of defining the property values is required by a 
specific problem, subroutines PKX , PKY , PROC or PQ which 
calculate respectively the conductivity in x and y 
directions, the thermal capacity and the heat generation per 
unit volume, can be easily modified or substituted. 

The shape functions used by this heat tranfer module 
are calculated by subroutines SHAPE and SHAP2. These 
routines are capable of constructing shape functions for a 
three-node triangle, a four-node linear quadrilateral, an 
eight-node quadratic serendipity quadrilateral, a nine-node 
quadratic Lagrangian quadrilateral or any combination 
i n-be tween. 

A four- to nine-node two dimensional element is 
shown in Figure 3. As shown by Bathe and Wilson in Ref. 5 
the temperature within the element is expressed at any time 
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in the local coordinate system s, t in terms of the nodal 
temperatures by 



T(s,t) = I N i (s / t) T i 



where 



N x = (l-s)(l-t)/4 - (N 5 +N g )/2 - Ng/4 

N 2 = ( 1+s ) ( 1-t )/4 - (N 5 +N 6 )/2 - Ng/4 

N 3 = ( 1 + s) ( 1+t )/4 - (N ? +N 6 )/2 - Ng/4 

N 4 = ( 1-s) (l+t)/4 - (N ? +Ng)/2 - Ng/4 

N s = ( 1-s 2 ) ( 1-t )/2 

N £ = (l+s)(l-t 2 )/2 
0 

N ? = (l-s 2 )(l+t)/2 
N g = ( 1-s) (l-t 2 )/2 
N g = (l-s 2 )(l-t 2 ) 

If any of the nodes from five to nine are omitted 
the corresponding value of N is zero. 

The connectivity table must follow the numbering 
convention shown in Figure 3, otherwise wrong results v/i 1 1 
be obtained. 

The maximum number of nodes per element will depend 
on the type of elements to be used and may be from three to 
nine. The eight-node serendipity shape functions occur if 
the ninth node number is omitted. The four-node 
quadrilateral shape functions are computed if only the first 
four nodal connections are non zero and the three-node 
triangle if the first three nodal connections are non zero. 
Final 1 y, if a mid-side nodal connection is omitted the edge 
is 1 i near. 
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The line boundary conditions, closely related to the 
type of element used, are defined in the cards following the 
macro MATE. The line subjected to a specified boundary 
condition is identified using the local coordinates. A line 
is numbered 1 or 2 if it is perpendicular to the axis s or 
t, respectively. Those numbers are positive or negative 
according to which direction of the axis it is 
perpendicular. Thus, as an example, the line defined by the 
nodes 2, 6 an 3 is numbered 1 while the line defined by 

nodes 1, 5 and 2 is numbered -2. 

The boundary conditions are treated in subroutine 
BC0ND2 and four types are considered: 

a) specified flux 

b) convection with constant heat tranfer coefficient 

c) convection with temperature dependent heat 
transfer coefficient 

d) radiation 

If more complex boundary conditions are required 
subroutine BC0ND2 may be easily modified. 

The numerical integration necessary to evaluate the 
matrices and vectors in equation (6) are performed using the 
Gauss quadrature formula; the abscissas and weight 
coefficients are calculated in subroutine PGAUSS. The user 
may choose the number of points per direction used in the 
integration from one to six, but the default value is four 
points . 

To illustrate the above, we assume that our sector 
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used as an example Is part of the cross section of a hollow 
cylinder made of a material with conductivity in x and y 
directions of 10, specific heat 250 and density 8000. If the 
outside surface is subjected to convection to an ambient 
atmosphere at 20 degrees with a constant heat tranfer 
coefficient of 200, the cards following the macro card MATE 
are 
MATE 

1 2 

10. 10. 8000. 250. 3 

621 200 . 20 . 

it must be noted that a plane analysis is specified 
and the integration will be performed using three points per 
direction. The reader must also note that this version of 
the problem Is slightly different from the one given when 
the mesh data preparation was discussed. In the previous 
problem the line now subjected to a convection boundary 
condition was at a specified constant temperature of 20. If 
this second problem i s to be solved the cards corresponding 
to nodes 31-33 in the BOUN and FORC sets must be omitted. 

Following is a list of subroutines included in the 
two dimensional heat tranfer module: 

ELMT02 : main routine; forms the necessary arrays 

for the solution. 

SHAPE : calculates the shape functions for the 

triangle and the linear quadrilateral elements. 

SHAP2 : adds the quadratic terms and the center node 
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to the shape functions 



PGAUSS : determine the abscissae and weight 

coefficients for the numerical integration. 

BC0ND2 : adds to the arrays calculated in ELMT02 the 
contr i bui tion from the specified boundary conditions. 

PKX : determine the temperature dependent 

conductivity in the x direction. 

PKY : determine the temperature dependent 

conductivity in the y direction. 

PROC : determine the temperature dependent thermal 

capac i ty. 

PQ : determine the temperature dependent heat 
generation per unit volume. 

CONV : determine the temperature dependent heat 

transfer coefficient. 

TABLE : called by PKX, PKY, PQ, PROC and CONV; 
calculates the correspondent coefficient to a given 
temperature by linear interpolation between two consecutive 
entr i es in a table. 

JACBB2 : calculates the jacobian determinant when 

an integration over a line is necessary. 

2 . Three Dimensional Heat Transfer Element 

This module is called ELMT03 and is accessed using 
the number three in column ten of each material number card. 

This element is the generalization of ELMT02 for 
three dimensional space and little remains to be said about 
its use. 
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The material properties are read in a similar way as 
the two dimensional case/ but the conductivity in the z 
direction must be added. This module uses the same 
subroutines as ELMT02 to read and calculate the temperature 
dependent properties. 

Subroutines SHAP3D and SHAP3 may construct 
interpolation functions for any combination between a 
eight-node linear brick and a 21-node Lagrangian brick. 

An eight- to twenty one-node three dimensional solid 
element is shown in Figure 4. The natural coordinates (r, s 
and t) of the eight corner nodes are (+1/ ±1/ ±1) / of the 
twelve mid-edge nodes are (0/ +1/ ±1)/ (±1/ 0/ +.1) and (+.1/ 
+1/ 0) and of the center node are (0/ 0/ 0). 

The temperature within the element is defined in 
terms of the nodal temperatures at any time by [Ref. 5] 

T(r/S/t) = Z N i (r /S/ t)T i 

If we define 



G( 0 


' e i> " 


.5(1+0^) / for 0.* +1 


and 


= 


l-B 2 for 6 i = 0 


g i 


= G( r# r 


^)G( S/ s dG( t, t^) 


where r^. 


s^ and 


t^ are the natural coordinates of 


element nodal points/ the interpolation functions are 


N i 


= g l 


- <N 9 .N 12 *N 17 )/2 - N 21 /8 


N 2 


S 2 


- <N 9 * N 10* N 18 )/2 - N 21 /8 


N 3 


s g 3 


- (N 10 *N 11 *M 19 )/2 - N 21 /8 


N 4 


= S 4 


- <N u .N 12 .N 20 )/2 - N 21 /8 
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n 5 ■ *5 * (N 13* N 16* N 17 )/2 * N 21 /S 

N 6 ' - <N 13* N 14* N 18 )/2 * N 21 /S 

N 7 *'«7 - (N 14* N 15* N 19 )/2 - N 21 /S 

N g - * 8 - !N 15 *N 16 *N 20 )/2 - N 21 /8 

Nj = gj - N 2 i/ 4 for j =9, . . . , 20 

N 21 = S 21 



If any of the nodes from nine to twenty one are 
omitted the corresponding value of g I s zero. 

The numbering convention used In the definition of 
the shape functions and shown in Figure 4 must be followed 
by the user when the connectivity table is built and the 
surface boundary conditions are specified. The surfaces are 
numbered .+1, +.2, +.3 as they are perpendicular to the 
positive or negative directions of the axis r, s, t, 

respecti vel y. 

As the numerical integration routine is the same as 
for ELMT02, the user may choose from one to six points per 
direction, the default value being four points. 

Following is a list of subroutines included in the 
three dimensional heat transfer module: 

ELMT03 : main routine; forms the matrices necessary 



for the solution. 



SHAP3D : calculates the shape functions for the 
eight node linear brick. 

SHAP3 : adds the quadratic terms to the shape 

functions and the center node for the Lagrangian brick. 

BC0ND3 : adds to the matrices the contr i bui tion from 
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the specified boundary conditions. 

PKZ : determine the conductivity in the z direction 
JACB63 : calculates the jacobian determinant when an 
integration over a surface is necessary. 

C. TIME INTEGRATION MODULE 

The first order ordinary differential equation solver is 
accessed when the macro command 0DE1 is utilized and uses 
the Zlenkiewicz two- and three-time level schemes. The 
details of these algorithms are treated in section IV. 

This module was designed for the solution of first order 
equations but can, with little programming effort/ be 
expanded in order to include higher order algorithms and 
solve higher order differential equations. 

In order to economize computer memory space at some cost 
of execution time, this module only uses the space 
previously reserved for the mesh construction. As the mesh 
data is not needed during the execution of a time step 
integration, all the corresponding arrays are kept in file 9 
during that period. Once the next displacement vector is 
calculated and the current displacement vector is saved in 
file 10, all the information in file 9 is retrieved and the 
mesh data restored. 

File 10 is used as a working space during the time 
integration and the current displacement vector is saved in 
it for possible future use. This vector is retrieved if the 
algorithm used In the next time step integration is a three 
point scheme. 
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Besides the time step size change using the appropriate 
macro instructions (macro TOL), an optional automatic time 
step adjustment was incorporated in the program. The norm of 
the difference between the displacements vectors at two 
consecutive times is computed at each step. If the norm is 
less than a predetermined value AT max f the time step size 
is doubled before going to the next step, whereas if the 
norm is greater than a prespecified value the time 
step size is halved and the calculation for that time step 
is repeated until the norm is acceptable. The magnitudes of 
the maximum and minimum values of the norm are problem 
dependent and must be choosen by the user. If they are not 
specified no time step adjustment will be performed. The 
recalculation of the displacement vector is allways done 
using the two point scheme, even when a three point scheme 
was utilized for the first calculation. 

The macro command 0DE1 used to access this module 
must be followed by a second macro command in order to 
determine the function to be performed. The secondary macro 
instructions are I N I T , LINE or QUAD. 

The couple 0DE1 INIT reads the integration constants 
theta , beta and gamma, the parameters for the automatic 
time step adjustment and the initial displacement vector. 
This data must follow the macro program (see user's 
instructions). No time integration is performed by this 
l nstruct l on. 

The couple 0DE1 LINE performs the two point scheme 
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and the current displacement vector is substituted by the 
new calculated displacement vector. 

The couple 0DE1 QUAD has a function similar to 0DE1 
LINE but uses the three point scheme. 

Once the execution of an 0DE1 macro command is 
complete/ the mesh data is restored and the displacement at 
time t replaced by the displacement at time t+ At/ when the 
secondary macro command is LINE or QUAD. If I NIT is used, 
the displacement vector becomes the given initial 
displacement vector. 

The subroutines included in this module are: 

P0DE1 : control routine; reserves working space and 
calls the appropriate subroutines. 

INIT : reads the constants theta, beta and gamma and 
the maximum and minimum values for the norm used in the time 
step adjustment. 

PLINE : performs the two point integration scheme. 

PQUAD : performs the three point integration scheme. 

TERM : performs the automatic time step adjustment. 
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IV. APPLICATION OF FEAP TO HEAT TRANFER PROBCEMS 



The program FEAP requires from the user the knowledge of 
the algorithm to be used In the solution of the problem to 
be treated. 

In this section the algorithms used In the solution of 
heat conduction problems are discussed. 

A. STEADY STATE PROBLEMS 

1 . Algorithms for Linear Problems with Macro Program 
In this case equation (8) becomes 

[K] {T} + {F} = {0} (12) 

The algorithm used to solve (12) Is described by 
[K] {v} = -{F} -[K] {T 0 } (13) 

{T} = {T° } + {v} (14) 

where {T 0 } are the specified nodal temperatures. 

The macro Instruction TANG builds the left hand side 
of (13) while FORM forms the vector In the right hand side. 
The Instruction SOLV solves the system (13) and performs the 
step defined by (14). 

Consequently the simplest macro program used to 
solve this problem Is 
TANG 
FORM 
SOLV 
Dl SP 
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2 . Algorithms for Nonlinear Problems with Macro 



Program 

Equation (12) still applies to this case but now [K] 
and {F} may be dependent on the nodal temperatures {T}. The 
well known Newton-Raphson iteration may be used and the 
algorithm is described as 

[KCT 1 )] {v 1 } = - {F (T 1 ) } - [KCT 1 )] {T 1 } = {R 1 } (15) 

{T i+1 } = { T 1 } + {v 1 } (16) 

where {T 1 } are the nodal temperatures at the ith iteration. 
The algorithm defined by (15) and (16) is the one used for 
the linear problem repeated several times. Then the macro 
program may be 
LOOP n 

TANG 
FORM 
SOLV 
Dl SP 
NEXT 
DISP 

where n are the user's guess of the number of iterations 
necessary to obtain equilibrium. However, the program has an 
internal check on the value of the vector norm | | R 1 1 I , 
where 

| |R | | = (Z |R j| 2 ) 1/2 

k=l K 

Whenever 

1 IR 1 1 1 < TOL x max | |R^ | | 

j 

where TOL Is a predefined tolerance ( 10 is the default 
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value), the iteration ceases and a skip to the macro command 
immediately following the first NEXT occurs. Usually one 
begins with zero as the initial guess of the displacement 
vector; however, any other vector may be used. 

B. TIME DEPENDENT PROBLEMS 

1 . Two and Three Point Recurrence Schemes for First 

Order Equation 

The set of differential equations 

[K] {T} + [ C] {T> + {F} = {0} (8) 

may be solved using one of the many recursive schemes 
described in the literature. From them, the two and three 
time level schemes presented by Zienkiewicz [Ref. 4] offer a 
wide range of choices for the solution of linear and 
nonlinear problems. 

The recurrence relation for the two-time level 
scheme can be wri tten as 

( ^[C] + 0[K]){T n+1 } + c-^rc] + (l-0)[K]){T n } + 

{F} = {0} , 0 < 0 < 1 (17) 

where the subscripts n and n+1 denote evaluation at time t 
and t+At. The vector {T} is defined as 

{?} = (F n+1 > + (l-0){F n > (18) 

The choice of the parameter 0 defines the particular 
scheme to be used and the reader will recognize a well known 
series of finite difference formulas with a modification of 
using a weighted loading term ("F} . 

Consider the system of decoupled equations in terms 
of the modal participation variables T^ 
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0 



(19) 



c. T. + k. T. + £. = 

11 li l 

for free response/ i.e./ f^O, expression (17) becomes 

<; h C i + k i e) < T iVl * c i * CT.J n -0 ( 20 ) 

Substituting the relation 



‘ X 'Vn (21) 

in (20)/ the resulting characteristic equation of the 
recurrent scheme solved for X gives 

X = [l-k^(l-0) At/c^] / [l+k^0 At/c^] 

The scheme will be unconditionally stable if 

I X | < 1 



for any value of p^ where 
V i = (kj/Cj) At 

This condition is satisfied for 

0 > 1/2 



( 22 ) 



On the other hand i f 

0 < 0 < 1/2 

stability is conditional requiring 

Pi < 2/ (1-20) 

Figure 5 shows how X varies with p^ for the schemes 
discussed later and how it compares with the exact value of 

X = exp(-p i ) 

The schemes considered are 

a) Crank-Ni col son with 0*1/2 

b) 0*2/3 proposed by Zienkiewicz [Ref. 4] 

c) 0=3/4 proposed by the author 

d) 0=0.873 proposed by Liniger in Ref. 6, 

The Zienkiewicz three-time level scheme is defined 
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as 



(Y [C] At[K]){T n+1 > + ( (1 - 2y) [C] + (l/2-2B+y) At [K] ) {T n > + 

(- (1-Y) [C] + (l/2+6-Y)At[K]){T n _ 1 } + At{F} = {0} (23) 

where the subscripts n, n+1 and n-1 denote evaluation at 
time t, t+At and t-At respectively. 

The vector { F} is defined as 

{F} = B{F n+1 > + (l/2-2B+y){F n } + (1/2 + B-y) {F^} (24) 

and the parameters y and g define the particular scheme to 
be used. 

Considering again the system of decoupled equations 
(19) and using the relations (21) and (22) one finds that 
the characteristic equation for the three-time level scheme 
i s 

(y+BPj) X 2 + [ (l-2y) + (l/2-26+y)P i ]X + 

[- (1-y) + (l/2 + B-y)p i ] = 0 (25) 

Wr i ting 

g = [l+(l/2+y)p i ]/[y+Bp i ] 
h = [-l+(l/2-y)p i ]/[y+Bp i ] 

the roots of (25) may be written as 

4 

X, ? = (2-g) / 2 ± [(2-g) 2 -4(l-h)] 1/2 /2 

1,2 

These roots will be complex if the quantity under 
the square root is negative. Then the modulus of X is 
I X | = (l-h) 1/2 

The scheme is stable if | X | <1 for any p^and Wood in 
Ref. 7 shows that this condition is satisfied if 

y>l/2 and B>y/2 
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Equation (25) has two roots: X 2 / the principal root 
which is the aproximation to the exact 

X = exp(-p i ) 

and the spurious root X 2 . For relative stability one must 
have 

U x l < U 2 I 

Also a negative root or complex roots can produce 
osci 1 lation. 

In Figure 6a-f the real roots X^ and X ^ and the 
modulus | X | for the complex roots are plotted against p^ 
for the schemes: 

a) 3*1/3/ y-1/2 proposed by Lees in Ref, 8 

b) 3=3/4/ y=1 proposed by Hogge in Ref. 9 

c) 3*0.646/ y=1.2184 proposed by Wood [Ref. 7] 

d) 3=4/5/ y = 3/2 proposed by Zlenklewicz [Ref. 4] 

e) 3=9/10/ y = 3/2 proposed by the author 

f) Fully implicit algorithm, 3=1, Y=3/2 [Ref. 7], 
From Figure 6a one may predict a very strong 

oscillatory behaviour for the Lees algorithm a). 

In order to study the various schemes, the physical 
problem represented by the nondimensi onal ized linear heat 
conduction equation 

5T = 6^T 
** 6x 2 

was solved in a bar of length 4 and width 1 subjected to the 
boundary conditions 



T=1 


at 


0 

n 

X 


f-0 

5x 


at 


x=4 
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and initial condition 

T=0 at t=0 

i.e., consider a step change in the surface temperature. 
This problem was also studied by Wood and Lewis in Ref. 10. 

The problem was discretized in space using ten equal 
length two dimensional quadratic elements through the length 
of the bar. This type of elements was chosen because it was 
reported by Wood and Lewis as giving the worst numerical 
resul ts . 

In order to avoid the discontinuity in the loading 
term caused by the step change of the surface temperature at 
the initial time, the problem was first solved starting from 
time t = At, where At is the time step size. The value of the 
temperature distribution at t=At is provided by the exact 
analytical solution. This procedure also provides the 
necessary two starting vectors for the three-time level 
schemes. 

The temperature of the nodes at x=l is used as 
reference value to compare the various schemes. 

in Figure 7a-d the results obtained from the 
various two-time level schemes with At=2 for the temperature 
at x=l are compared with the analytical solution. For this 
Ideal starting conditions the Crank-N i col son, 0=1/2, proves 
to be the most accurate algorithm. When the step size was 
Increased to 10, all the schemes performed well as shown in 
Figure 3a-d. 

The corresponding results for the three-time level 
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schemes are presented in Figure 9a-f and Figure lOa-f for 
At=2 and At=10, respectively. 

Strong oscilations were observed with the Lees 
algor i thm a) . 

For At=2 the schemes c), d) and e) performed well. 
For the larger step size of 10 the safest schemes are e) and 
f). 

The step change in the loading term has been 
reported [Ref. 10] as producing relative instability when 
some of the schemes are used, particularly with the 
Crank-Ni col son algorithm. The noise may be supressed 
reducing the magnitude of the time step but this is 
impraticable in most problems. Other techniques may be used 
to reduce the amplitude of the oscillations as those 
discussed by Wood and Lev/i s [Ref. 10] and Gresho and Lee in 
Ref. 11. 

The step change in the loading term in the problem 
in study is due to the variation of the surface temperature 
at the initial time t=0. That discontinuity In the force 
term Is Illustrated In Figure 11a, When the time dimension 
is discretized it has been common practice to calculate the 
numerical solution at the nodes t=0, t=At, etc. When a 



three-time 


level 


scheme i s 


used 


the node at t=- At 


i s 


cons i dered 


assumi ng 


al 1 


condi t i ons 


steady before t=0. 


The 


force values used 


by 


this 


method 


are indicated by 


the 


cl rcl es in 


Figure 


11a, 


This 


procedure transfers to 


the 



numerical solution the uncertainty in the value of the force 
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term at time t®0 introduced by the mathematical 
discontinuity assumed In the analytical solution. This 
problem is avoided if one uses the procedure illustrated in 
Figure lib. The nodes considered for the discretization are 
those at t=-At/2, t=At/2, t=3At/2 / etc. The discontinuity at 
t=0 is smoothed by the interpolation of the loading term 
inherent to the particular scheme being used. The values of 
the force term for every node are well defined and no 
uncertainty is associated with any of them. For the 
three-time level schemes the node at t=-3At/2 is used 
assuming all conditions steady before t*0. 

The problem studied was solved using both starting 
procedures. The results obtained wl th the first procedure, 
start at t=0, using the two- and three-time level schemes 
are shown in Figures 12a-d and 13a-f, respectively. The 
second procedure, start at t=- At/2, gives the results shown 
In Figures 14a-d and 15a-f for the two- and three-time level 
schemes, respectively. 

One may conclude that, in this simple problem, when 
the procedure illustrated in Figure lib is followed, the 
step change in the loading term do not deteriorate the 
accuracy of the numerical results obtained with the schemes 
tested, with the exception of the Crank-Nicol son algorithm. 
This scheme shows an oscillatory behaviour not observed when 
the step change In the force term Is not present. 

The starting procedure of Figure lib also proves to 
be an efficient way of starting the time integration with 
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the three-time level schemes. 

2 . Algorithms for Time Dependent Problems with Macro 
Program 

A simplified form of the two-time level scheme is 
available in FEAP and is described by 

* e [ K n1>'W + * Cl*e) [KJ 

At{F} = {0} (26) 

where the subscripts n and n+1 denote evaluation at time t 
and t+At, respectively. 

The three-time level scheme programmed in FEAP is 
described by the expression 

( Y [C n ] +6 At[K n ]){T n+1 ) + ((l-2y)[C n ] + (1/2-20+y) At [KJ ) {T n > 
+ C-(l-Y)[C n ] + (l/2 + 6-y)At[K n ]){T n _ 1 } 

+ {F n > = {0} (27) 

where the subscripts n, n+1 and n-1 denote evaluation at 

times t,At + t and t-At respectively. 

It must be noted that in (26) and (27) the vector 

{F }is an aproximation of the Interpolated force vector {F}. 
n 

If the force vector is strongly dependent on time, 
this aproximation may produce wrong results. 

In problems where the stiffness matrix is constant 
and the force term results from the specified boundary 
displacements or where the specified boundary forces are 
only time dependent, the interpolation of the force term may 
be easily acompllshed since the force/displacement boudary 
values can be changed at any time using the macro commands 
MESH or PROP. This procedure was found partcularly useful 
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when a step change In boundary displacements or forces Is 
applied at the initial time. 



The first order ordinary differential equation 
solver module is acessed using the macro command 0DE1 
followed by one of the following macro commands: 

INIT to specify the Initial displacement vector 
and the integration constants 
LINE to perform the two-time level algorithm 
QUAD to perform the three-time level algorithm 
Since the matrices [K] and [C] and the vector {F} 
must be formed before the time integration/ the macro 
commands TANG/ CMAS or LMAS and FORM are closely associated 
with the use of 0DE1. It should be noted that the macro FORM 
forms the vector 

< R<T n>> ■ - {F n> - t K nl <V 

which is allways dependent on the current displacement. 
Therefore the macro command FORM must be used every time 
step and preceed 0DE1, 

To solve a fully nonlinear problem the following 
macro program may be used: 

DT At 

0DE1 INIT 
LOOP n 

TANG 

CMAS ( or LMAS ) 

FORM 

0DE1 LINE ( or QUAD ) 
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TIME 



DISP 

NEXT 

The macro command program must be followed by the 
data cards necessary to specify the integration constants 
6/ y and 3 and the initial displacement vector. 

If the matrices [K] and/or [C] are not temperature 
or time dependent then the instructions TANG and/or CMAS or 
LMAS must be placed outside the loop. 

For the simplest case of a linear problem, the macro 
program may be 
DT At 

0DE1 I N IT 
TANG 

CMAS ( or LMAS ) 

LOOP n 

FORM 

0DE1 LINE ( or QUAD ) 

TIME 

DISP 

NEXT 



40 



V. NUMERICAL EXAMPLES 



A. RADIAL TEMPERATURE DISTRIBUTION IN A HOLLOW 

CYLINDER 

Consider the hollow cylinder of infinite length already 
discussed in section III. Consider the case in which the 
inside wall is maintained at a constant temperature and heat 
is lost by convection through the outer wall. This problem 
was considered for solution by Lew [Ref. 3]. 

The geometrical and thermophysical character i sties are: 



Inside radius ( r ^ n ) 


0.1 


m 


Outside radius ( r ou t) 


0.3 


m 


Thermal conductivity (k) 


10 


W/m °C 


Outside ambient temperature 


20 


°C 


Heat transfer coefficient 


200 


W/m 2 °C 


Inside wall temperature (T^ ) 


820 


°C 


The finite element model of this problem 


was completly 



discussed before and is shown in Figure 1. The complete 
computer output is given in Appendix B. It must be noted 
that the element arc width used is arbitrary since the 
problem is symmetric and no heat is conducted perpendicular 
to the radius. 

The comparision between the analytical and FEAP 
solutions is shown in Figure 15. In this case the values 
provided by the finite element model used are very accurate 
and no further mesh refinement is necessary. 
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B. NONLINEAR STEADY STATE HEAT CONDUCTION IN A SLAB 
WITH TEMPERATURE DEPENDENT CONDUCTIVITY 

A liquid is boiled by a flat electric heater plate of 

thickness 2L. The internal energy Q generated electrically 

may be assumed to be uniform. The boiling temperature of the 

liquid/ corresponding to a specific pressure, is T w . 

In the finite element analysis for the temperature 

di str i bu i ti on in the heater, ten equal length quadratic 

elements were used to represent the unit cross section 

through the half plate thickness L, since the problem is 

symmetric. The extremity of the resulting slab is assumed at 

the temperature T while the other is insulated. 

w 

The thermal conductivity is assumed temperature 
dependent according to the expression 

k = a( 1+bT) 

where a and b are constants. 

Assuming the following values for the parameters in a 
consistent system of units 

T =0.0 

w 

a = 0.5 
L = 1.0 
Q = 1.0 

the problem was solved for b equal 0.0, 1.0 and 2.0. 

The numerical results are compared with the analytical 
solutions in Figure 16 and one may conclude that FEAP 
provides exact solutions in this problem. 
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C. TRANSIENT TEMPERATURE D I STR I BU IT I ON IN A HOLLOW 
CYLINDER 

The same hollow cylinder of infinite length treated 

before is considered again but the case solved is the one 

with both outside and inside walls maintained at constant 

temperatures. Thermal diffusivity/ a=k/ pc is specified as 
2 

0.00005 m /sec. Initially the cylinder is at an ambient 
temperature Tq of 20°C; suddenly the temperature of the 
inside wal 1 is raised to 820°C while the outside wall is 
maintained at 20°C. All other conditions are the same as for 
the steady state problem. 

In order to obtain the unsteady solution the number of 
elements in the finite element model was doubled# thus 
twelve radial quadratic elements were used. The element arc 
width chosen is the same 6° as before. 

The two time-level scheme with 9=2/3 was used for the 
time integration. In order to compare the results the 
dimensionless temperature T*# radius r* and Fourier number 
Fo were used. They are defined as 

T* - (T - T 0 )/(T ln - T q ) 

r * * r/r in 

and 

Fo = at/r? n 

A time step size At of 2 seconds was chosen initially. 
After 10 sec# At was increased to 10 sec and to 50 sec after 
200sec. After 1000 sec#At was increased to 100 sec. 

Figure 17 shows the evolution of the temperature of the 
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center interior points of the cylinder and Figure 18 
compares the analytical values with the FEAP solution for 
the radial temperature di str i bu i t ion for two different 
times. 

D. SLAB WITH RADIATION-CONVECTION BOUNDARY CONDITIONS 

Consider a plane slab of thickness L, with constant 
thermophysical properties (k= p 0 = 1 )/ subjected to 
simultaneous radiation and convection at x=0 and insulated 
at x = L (0<x<L). The slab is initially at a uniform 
temperature T Q and it was decided to consider zero ambient 
temperatures for convection and radiation. 

It was also decided to consider the Biot number 
Bi = h c L/k = 1 

where h is the convection heat tranfer coefficient and 
c 

R = eaTpL/k = 4 

where e is the emissivity and a the Stefan-Bol tzman 
constant. 

The same mesh as in Example B was employed for the 
finite element solution. The time integration was performed 
by the three-time level scheme with y =1/2 and 6=1/3. A 
variable time step size was chosen. A value of At=0.001 was 
used initially. After t=0.02/ At was increased to 0.0025. 
After t=0,2/ At was increased to 0.01. After t=0.4/ At was 
increased to 0.025. After t=1.0/ At was increased to 0.05. 

The history of the ratios T^/T^ and T^/Tq are plotted in 
Figure 19 against the Fourier number 
Fo = at/L^ 
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T is the surface tmperature at x = 0 and T is the surface 
w m 

temperature at x=L. 

The finite element results are compared with those given 
by Haji-Sheick and Sparrow in Ref. 12 and obtained from a 
Monte Carlo method. 
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VI. CONCLUSIONS AND RECOMMENDATIONS 



The computer program in this thesis provides an accurate 
and reliable means for solving a variety of Heat Transfer 
problems. The use of this program and efforts to increase 
its versatility are highly encouraged. 

The capabilities of the program, while quite 
significant, can still be improved. The present version is 
designed for an "in-core" solution technique, which 
restricts the problem size to within the computer core-size. 
The capacity to handle large problems may be increased 
through the use of some "out of core" technique for solving 
the system of equations and even for building the mesh data. 
In that case the size of the problems treated would be 
restricted only by the availability of external storage 
devices. 

In order to check the mesh input and to easily interpret 
the output, a graphics module must be included in the 
program. 

Although the time integration module seems to provide 
reliable results in linear and some nonlinear problems, it 
should be subjected to further research in order to test it 
with different types of nonlinear analysis. 
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Figure 1. Mesh for Hollow Cylinder Problem 




Figure 2. Eight -Node Serendipity Quadrilateral 
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Figure 3. Local Node Number Sequence for a Quadratic 
Lagrangian Element 
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Figure 4. Local Node Number Sequence for a 21 -Node 
Three Dimensional Element 
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Figure 5. Zieiikiewicz Two-Time Level Scheme 
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Figure 6 . (continued) 
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Figure 7. 
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Temperature Temperature 



l.l 




(c) 




(d) 

Figure 7. (continued) 
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Temperature Temperature 
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Figure 8. 
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Temperature Temperature 




(c) 




(d) 

Figure 8 . (continued) 
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Temperature Temperature 
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Cb) 

Figure 9. 
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Figure 9 . (continued) 
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Temperature Temperature 
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(e) 




(f) 

Figure 9. (continued) 
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Temperature Temperature 




(a) 




(b) 



Figure 10. 
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Temperature Temperature 





(d) 

Figure 10 . (continued) 
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Temperature Temperature 




(e) 




(£) 

Figure 10 . (continued) 
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(a) 




Figure 11. Step Change in the Loading Term 
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Temperature Temperature 
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Cb) 



Figure 12. 
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Temperature Temperature 




(c) 




Cd) 

Figure 12. (continued) 
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Temperature Temperature 




(a) 




Cb) 

Figure 13. 



67 



Temperature Temperature 
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Figure 13. (continued) 
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Temperature Temperature 
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Figure 13. C con tinued) 
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Temperature Temperature 




(a) 




(b) 

Figure 14. 
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Temperature Temperature 




(c) 




(d) 



Figure 14. (continued) 
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Temperature Temperature 




(a) 




Cb) 

Figure 15. 
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Temperature Temperature 
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(c) 




(d) 



Figure 15. (continued) 
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Temperature Temperature 



l.l 





(f) 

Figure 15. (continued) 
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Steady State 
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Figure 16. Radial Temperature Distribution in a Hollow Cylinder. T 
temperature T and the radius r are expressed in degrees 
Celsius and meters, respectevely . 
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Figure 17. Effect of Linear Conductivity Variation on the Temperature 
of a Plate with Uniform Heat Generation. 
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Figure 18. Unsteady Temperature of the Center Interior Points in 
Hollow Cylinder. 
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Figure 19. Unsteady Radial Temperature Distribution in a Hollow Cylinder 
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Figure 20. Unsteady Temperature Distribution in a Slab with 
Boundary Conditions. 
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APPENDIX A 

USER INSTRUCTIONS FOR FEAP 
A • 1 TITLE ANC CONTROL INFORMATION 



TITLE CAFC-FCRMAT(2QA4) 

Tf-5 TITLE CARC ALSO SERVES AS A START CF PROBLEM 
THE FIRST FOUR (4) COLUMNS MUST CONTAIN THE STAR 
fcCRD FEAF 

COLUMNS DESCRIPTION VARIABLE 

1 TC 4 MUST CONTAIN FEAP TITL(l) 

5 TO EC ALPHANUMERIC INFORMA- TITL(I),I 
TION TO EE PRINTED 
WITH OUT FUT AS PAGE 
HEADER. 



.... CONTROL CARC- FORMAT ( 71 f } 



COLUMN 


c 


DESCRIPTION 


VAPIAELE 


1 


TC 


G 


NUMBER OF NCCES 


NUMNP 


6 


TC 


1C 


NUMBER CF ELEMENTS 


NUMEL 


11 


TG 


15 


NUMBER CF MATERIAL SETS 


NUMMAT 


16 


TC 


20 


SPATIAL CIMENSION 
(MAXIMUM, UP TO 3) 


NDM 


21 


TC 


25 


NUMBER GF UNKNOWNS 
FER NOCE (MAXIMUM, 
UP TO 6) 


NDF 


26 


TC 


20 


NUMBER CF NOCES PER 
ELEMENT (MAXIMUM) 


NEN 


21 


TC 


■; e 


NOT USEC 





CARC. 



= 2,20 
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A. 2 CATA INPUT : MACRO CONTROL STATEMENTS 

INPUT RAC FC CONTROL CAPCS-F0RMAT(A4) 

THE INFLT CF EACH DAT* SEGMENT IS CONTROLLED EY THE 
VALUE ASSIGNED TO CC. THE FOLLOWING VALUES ARE ADMI- 
SSIBLE AND EACH CC CAPC MUST BE IMMEDIATELY F GLLOWEC 
EY THE APPROPRIATE DATA. 



CC VALLE 

COOR 

FCLA 

SPHE 

ELEM 

RATE 

EOLN 

FORC 



TERP 

PRIN 

NOPR 

PAGE 

ENC 



DATA TC eE INPUT 
COORDINATE DATA 
CONVERT PCLAP TO CARTESIAN 
COORDINATES 

CONVERT SPHERICAL TO CARTESIAN 

COORDINATES 

ELEMENT DATA 

MATERIAL CATA 

NODAL ECUNCARY CONDITION CATA 
PRESCRIBED NODAL DISPLACEMENT/FORCE 
DATA. FOR HEAT TRANSFER PROBLEMS 
THIS IS PRESCRIBED NODAL TEMFERA- 
TURE/FLL'X CATA 
TEMPEF ATL’PE DATA 
PRINT SLBSEQUENT MESH DATA 
(DEFAULT MODE) 

DO NCT PRINT SUBSEQUENT MESH CATA 
PRINTED OUTPUT CONTROL 
MUST eE LAST CARO IN MESH DATA, 
TERMINATES MESH INPUT 



EXCEPT FCF THE ENC CAPC THE DATA SEGMENTS CAN BE 
IN ANY CPDER. IF THE VALLES OF BOUN,FGRCt CF TEMP 
ARE ZEFC » NC INPUT CATA IS REQUIRED. 



COCRDI NATE CATA - FORM AT (2 15 , 7F10 .0 ) 

MUST IMMEDIATELY FOLLOW A COOR MACRO CARD. 

THE CCCFCIN ATE CATA CARD CONTAINS THE NODS NUMBER N 
AND THE VALLE OF THE COORDINATES FOR THE NCCE. CNLY 
THE VALUES CF (XL ( I ) , 1=1 , NCM ) ARE USED, WHERE NDM IS 
THE VALLE INFLT ON THE CONTROL CARD. 

NOCAL COGP C INATES CAN EE GENERATED ALONG A STRAIGHT 
LINE DESCFI EEC EY THE VALUES ON TWO SUCCESSIVE 
CARDS. THE VALUE OF THE NCCE NUMBER IS CCMPUTEC 
USING THE N ANC NG ON THE FIRST CARD TO COMPLLTE 
THE SEQUENCE N, N+NG,N+2*NG, ETC. NG MAY BE INPUT 
AS A NEGATIVE NUMBER, IF IT HAS INCORRECT SIGN THE 
SIGN WILL EE CHANGED. NOCES NEED NOT BE IN ORCER. 



COLUMNS DESCRIPTION 

1 TC 5 NCDE NUMEER 

6 TC 10 GENERAT CF INCREMENT 

11 TC 2C XI COORDINATE 

21 TC 3C X2 COOPCINATE 

31 TC <0 X3 COORDINATE 

TERMINATE WITH eLANK CARD(S) 



VARIABLE 

N 

NG 

XL( 1) =X(1,N) 
XL(2l=X(2,N ) 
XL ( 3 1 =X(3,N) 



.... POLAR TC CARTESIAN CONVERSION C ATA-FORMAT (3 15 , 2F10 . 0 ) 
MUST IMMEDIATELY FOLLOW A POLA MACRO CARD. 

THE FCLAF CCCRCINATES (R, THETA) MUST BE PREVIOUSLY IN- 
PUT AS CCCRC INATE DATA (MACRC CARD COOP) WHERE R IS XI 
ANC THETA IS X2. THE NOCAL POLAR COORDINATES CF NODES 
NI, NI+INC, NI+2*INC,..., NE ARE CONVERTED TC 
CARTESIAN COORDINATES. 

COLUMNS DESCRIPTION VARIABLE 

1 TC 5 FIRST NCCE TC BE NI 

CONVERTED 

6 TC 10 LAST NODE TO BE NE 
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11 


TC 


15 


CCNVERTEC 
NOOAL INCREMENT 




INC 


16 


TC 


25 


X COORDINATE 


OF 


THE 


XO 


26 


TC 


3 c 


CRIGIN 

Y COORDINATE 


OF 


THE 


YO 




CRIGIN 

TERMINATE WITH ELANK 


CARC ( S) 





.... SPHERICAL TC CART. CONVERSION DATA-FORMAT< 31 5 , 2F10 . 0 ) 
MLST I MMECC IATELY FCLLCW A SPHE MACRO CARC. 

THE SPHERICAL COORDINATES (R.THETA.PHI ) MUST BE PRE- 
VIOUSLY INPUT AS COORDINATE DATA (X1.X2.X2) FOLLOWING 
MACRC CAP C CCCP. THETA IS THE ANGLE FORMED BY R ANC X 
AND PHI IS THE ANGLE FCRMEO BY R AND Z. 

THE NCCAL SPHERICAL CCCPCINATES OF NOOES NI. NI+INC. 
Nl42*INC,....NE ARE CONVERTED TO CARTE SI A N CCCPDINATES 



COLUMNS DESCRI FTICN VARIABLE 

1 TC 5 FIRST NCCE TC BE NI 

6 TC 10 LAST NCCE TC BE NE 

CCNVERTEC 

11 TC 15 NODAL INCREMENT INC 

16 TC 2C X COORDINATE OF THE XO 

ORIGIN 

26 TC 25 Y COORDINATE OF THE YO 

ORIGIN 

26 TC 45 Z COORDINATE OF THE ZO 

CRIGIN 



TERMINATE WITH ELANK CARC(S). 



.... ELEMENT CATA - FORMAT(16I5) 

ML ST IMMEDIATELY FOLLOW AN ELEM MACRO CARC. 

THE ELEMENT CATA CARD CONTAINS THE ELEMENT NUMBER » 
MATERIAL SET NUMBER (WHICH ALSO SELECTS THE ELEMENT 
TYPE. SEE MATERIAL PROPERTY DATA) . AND THE SEQUENCE 
CF NCCES CONNECTED TO THE ELEMENT. IF THERE ARE LESS 
THAN NEN NCCES EITHER LEAVE THE APPROPRIATE FIELDS 
BLANK CR PUNCH ZEROS. 

ELEMENTS MUST BE IN ORDER. IF ELEMENT CARDS ARE 
CMITTEC THE ELEMENT 0*TA WILL BE GENERATEC FROM 
THE PREVIOUS ELEMENT WITH THE SAME MATERIAL NUMBER 
ANC ThE NCCES ALL INCREMENTED BY LX ON THE PREVIOUS 
ELEMENT. GENERATION TC THE MAXIMUM ELEMENT NUMEER 
CCCUPS WHEN A BLANK CARC IS ENCOUNTERED. 



COLUMNS 
1 TC £ 

6 TC 10 
11 TC 15 
16 TC 20 
ETC. 
ETC. 
ETC. 



DESCRIPTION 
ELEMENT NUMBER 
MATERIAL SET NUMBER 
NOOE 1 NUMBEF 
NODE 2 NUMBEF 

NODE NEN NUMEER 
GENERATION INCREMENT 



TEFMINATE WITH ELANK CARO(S). 



VARIABLE 

L 

IX(NENl.L) 
I X ( 1 . L ) 

I X < 2 * L 1 

IX (MEN . L ) 
LX 



.... MATERIAL PROPERTY DATA - FORMAT ( 15 .4X , II . 17A4 ) 
MUST IMMEDIATELY FOLLOW A MATE MACRO CARD. 

COLUMNS DESCRIPTION VARIABLE 

1 TC 5 MATERIAL SET NUMBER MA 

6 TC S 

10 ELEMENT TYPE NUMBER I EL 

11 TC TE ALPHANUMERIC INFORMATION XHEC 

TO BE OUTPUT 

EACH MATERIAL CARD MUST BE FOLLOWED IMMEDIATELY 
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BY THE MATERIAL PROPERTY DATA REQUIRED FCR THE 
ELEMENT TYPE I EL BEING USED, E. G. , SEE SECTION A. 4 
ANC A. 5 FCR HEAT TRANSFER ELEMENTS. 



ECUNCARY RESTRAINT DATA - FORMAT ( 16T 5) 

MUST I MME D IATELY FOLLOW A EOUN MACRO CARD. 

FOR EACH NCDE WHICH HAS AT LEAST ONE DEGREE CF FREEDOM 
WITH A SPECIFIED DISPLACEMENT, A BOUNDARY CONDITION 
CARD NIST EE INPUT. THE CONVENTION USEC FCR BOUNDARY 
RESTRAINTS IS 

•EC .0 NO RESTRAINT, FORCE SPECIFIED 
.NE.C RESTRAINED. DISPLACEMENT SPECIFIED 
VALUES CF FORCE OR DISPLACEMENT INPUT IN FORC 



COLUMNS 
1 TO 5 
6 TO 10 
11 TO 15 

'ME . 20 

ETC. 



DESCRIPTION 
NODE NUMBER 
GENERATION INCREMENT 
DOF 1 BOUNDARY CODE 
DOF 2 ECUNCARY CODE 

COF NDF ECUNCARY CODE 



VARIABLE 

N 

NX 

I DL (1 I =ID(1 , N) 
I DL ( 2 ) - ID( 2 » N ) 



WHEN GENERATING 
SUESECUENT NODES 
INPUT .C-E.G, AND 
TIVE. ALL DEGREE 
CODES ARE ASSUME 
TERMINATE WI 



I DL ( NDF ) = 

ID ( NDF , N ) 
BOUNCARY CONDITION CODES FCR 
ICL IS SET TO ZERO IF IT WAS 
IS SET TO -1 IF INPUT NEGA- 
S OF FREEOOM WITH NON-ZERC 
D FIXEC. 

TH ELANK CARD(S). 



NOCAL FORCE C BOUNDARY VALUE DATA - FORMA T ( 21 5 ,7F10 . 0 ) 
MUST I N ME 0 IA TELY FOLLOW A FORC MACRO CARC. 

FOR EACH NCDE WHICH HAS A NON- ZERO NODAL FORCE CR 
DISPLACEMENT A FORCE CARC MUST BE INPUT OR GENERATED. 
GENERATION IS THE SAME AS FCR CO-ORDINATE CATA. THE 
VALUE SPECIFIED IS A FORCE IF THE CORRESPODINC- RES- 
TRAINT CODE IS ZERO ANC A DISPLACEMENT IF THE CORRES- 
PONDING RESTRAINT CODE IS NON-ZERO. 



COLUMNS DESCRIPTION 

1 TO 5 NODE NUMBER 

6 TC 10 GENERATION INCREMENT 

11 TO 2C DCF 1 FORCE (CISP. ) 

21 TO 20 CCF 2 FORCE (CISP.) 

ETC. 

ETC. DCF NDF FCFCE (DISP.) 

TERMINATE WITH A ELANK CARD. 



VARIAELE 

N 

NG 

XL(1)=F(1,N) 

XL(2)=F(2»N) 

XL (NDF ) = 

F ( NDF ,N) 



TEMPERATURE CATA CARD - FORMAT ( 21 5,F 10. 0 ) 

MUST IMMEDIATELY FOLLCW A TEMP MACRO CARC. 

NOT USEC IN HEAT TRANSFER PROBLEMS. 

FCR EACH NCDE WHICH HAS A NON-ZERO TEMPERATURE THE 
VALUE MUST EE INPUT. GENERATION OF VALUES CAN BE 
FERFCRMEC AS DESCRIBEC FOR CO-ORDINATES. 



COLUMNS DESCRIPTION 

1 TC 5 NODE NUMBER 

6 TC 10 GENERATION INCREMENT 

11 TC 2 C NODAL TEMPERATURE 

TERMINATE WITH ELANK CARD(S), 



VARIAELE 

N 

NG 

X L ( 1 I =T ( N ) 



FAGE CCNTRCL DATA - FCFMAT ( A 1 ) 

MUST IMMEDIATELY FOLLOW A PAGE MACRO CARD 
THE VALUE IN CCLUMN 1 CCNTFOLS THE PRINTED 



OUTPUT 
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ACCORDING TC Tl-E FOLLOWING CONVENTION: 

1 TITLE IS WRITEN AT THE TOP OF A NEW PAGE 
0 THE PRINTED OUTPLT IS CONTINUOUS WITHOUT 
SKIPPING PAGES ( DEFAULT MODS). 
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A. 2 PROBLEM SOLUTION: MACRC PROGRAMMING CCMANCS 



MACRO PROGRAMMING COMMANDS - F0RMAT< 2 ( A4 , IX) , F10.0 ) 
FOLLOWING IS A LIST OF MACFO INSTRUCTIONS WHICH 
CAN BE LSEC TO CONSTRICT SCLUTION ALGORITHMS. THE 
FIRST INSTRLCTION MUST BE A CARD WITH MACR IN 
COLUMN S 1 TC 4. THE INDICATED ISW VALUE IS USED BY 
EACH ELEMENT ROUTINE TC FEFFCRM THE APPROPRIATE 
CCMPUT AT IONS • 



COLUMNS 
1 TC 4 
CHEC 
CM AS 

CCN V 
DATA 



CISF 

**CT 

EIGE 



COLUMNS COLUMNS 
6 TC 9 11 TO 20 



EXCD 



FGPM 

LMAS 

LCCF 



MESH 



NEXT 

CCE1 



0CE1 

CCE1 

PRCP 

RE AC 
SOLV 

STPE 



TANG 

TIME 

*TIM 

**TCL 

UTAN 

ENC 



N 



N 



IN IT 



LINE 

CUAC 



N 



DESCRIPTION 

FERFCRM CHECK OF MESH (ISW=2) 
CONSISTENT MASS FORMULATION 
(ISW=5) 

DISPLACEMENT CONVERGENCE TEST 
READ DATA **MACRC COMMANDS. A 
**MACRC CARD IS INSEPTEC AS CATA 
FOLLOWING THE MACRC PROGRAM 
PRINT NOCAL DISPLACEMENTS EVERY 
N STEPS IN LOOP 
SET TIME INCREMENT TC VALUE V 
COMPUTE DOMINANT EIGENVALUE AND 
VECTOR OF CURRENT MASS AND 
SYMMETRIC STIFFNESS 
EXPLICIT CENTEREC DIFFERENCE 
INTEGRATION OF EQUATIONS CF 
MCTICN USING LUMPEC MASS. FIRST 
CALL RESERVES MEMORY ONLY. 

FORM RIGHT HAND SIDE OF EQUA- 
TIONS ( I SW=6 ) 

LUMPED MASS FORMULATION (ISW=5) 
LOOP N TIMES ALL INSTRUCTIONS 
BETWEEN MATCHING NEXT 
INSTRUCTION. 

INPUT MESH CHANGES (MUST NOT 
CHANGE BOUNDARY RESTRAINTS) .CA T A 
FOLLOWS MACRO PROGRAM. SEE 
SECTION A. 2. 

ENC CF LOOP INSTRUCTION 
INPUT INITIAL CONDITIONS AND 
CONSTANTS OF INTEGRATION FOR 
FIRST ORDER ORDINARY DIFFEREN- 
TIAL EQUATION SOLUTION (CATA 
FOLLOWS MACRO PROGRAM ) 

TWO-FOINT INTEGRATION SCHEME 
THREE-POINT INTEGRATION SCHEME 
INPUT PROPORTIONAL LOAD TABUS 
(CATA FOLLOWS MACRO PROGRAM) 
COMPUTE NOCAL RE ACT ICNS ( ISW = 6 ) 
SOLVE TANGENT EQUATIONS. UPDATE 
NOCAL DISPLACEMENTS. 

FPINT ELEMENT VARIABLES (E.G., 
STRESSES) EVERY N STEPS IN LOOP 
( ISW=4) 

SYMMETRIC TANGENT STIFFNESS 
FORMULATION (ISW=3) 

ACVANCE TIME BY CT VALUE 
INDICATE EXECUTION TIME FOR 
EVERY MACRO INSTRUCTION 
SET SOLUTION CONVERGENCE TOLE- 
RANCE TO VALUE V 
UNSYNMETRIC TANGENT STIFFNESS 
FORMULATION (ISW=3) 

ENC OF MACRO PROGRAM COMMANDS. 
CATA FOR PROGRAM FOLLOWS IN 
CRDER OF USE 



.... FIRST CRCEF C.C.E. SOLVER CATA 
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MUST FCLLCVi THE MACRO CARD END AND CORRESPONDS TO 
THE MAC PC I N STRUCT I CN CDE1 INIT. THE FIRST CARD 
INPUTS THE CONTROL INFORMATION AND THE NODAL INITIAL 
DISPLACEMENTS ARE INPLT ON THE FOLLOWING CARDS. 

NO AUTOMATIC TIME STEF AC JLSTMENT IS PERFORMEC IF 
CUMAX CR CLMIN ARE LEFT BLANK OR ZERO. 

.... CONTROL CARD - FORMAT (4F10.0) 



COLUMN 


c 


DESCRIPTION 


VARIABLE 


1 


TO 


1C 


INTEGRATION PARAMETER 
THETA FCP TWO PCINT 
SCHEME (DEFALLT VALUE 
THETA=2/2 ) 


C 5 


11 


TC 


2C 


INTEGRATION PARAMETER 
GAMMA FOR THREE POINT 
SCHEME 

(DEFAULT GAMMA=1 . 5 ) 


Cl 


21 


TC 


•C 


INTEGRATION PARAMETER 
BETA FCR THREE PCINT 
SCHEME 

(DEFAULT EETA=.8 ) 


C 2 


lit 

I - 1 


TO 


40 


MAXIMUM DISPLACEMENT 
ALLOWED 


DUMAX 


41 


TC 


50 


MINIMUM DISPLACEMENT 
ALLOWED 


DUMIN 



INITIAL CCNCITION DATA - FORMAT ( 2I5.7F10 .0 ) 

FOR EACH NODE WHICH HAS A NON-ZERO INITIAL DISPLACE- 
MENT A CAPO MUST BE INPUT CR GENERATED. THE GENERATION 
IS THE SAME AS FOR CCCFCINATE CATA. 



COLUMNS DESCRIPTION 

1 TO 5 NODE NUMBER 

6 TO 1C GENERATION INCREMENT 

11 TO 10 CCF 1 DISPLACEMENT 

ETC. 

ETC. COF NDF DISPLACEMENT 



VARIABLE 

N 

NG 

XL(1)=U(1+N-1) 
XL ( NCF ) = 



TERMINATE WITH BLANK CARD(S). 



.... FRCP CRT 1 0 N A L LOAD CARC - FORMAT C2I5.6F10 .0 ) 

A SIMPLIFIED PROPORTIONAL LOADING IS PERMITTED WITH 
PPOF = A1 + A2+TIME 4 A2*( S IN ( A4*T I ME 4 A5))**L 
WHERE THE COEFFICIENTS APE INPUT ON A CATA CARC 
FOLLOWING THE END MACRC CARC ACCORDING TO THE 
FOLLOWING TAELE. 

COLUMN DESCRIPTION 



6 


TC 


1C 


L 










11 


TC 


2 C 


MINIMUM TIME 
COMPUTEC 


FOR 


WHICH 


PRCP 


IS 


21 


TC 


(tl 

r% 


MAXIMUM TIME 
COMPUTED 


FCR 


WHICH 


FRCP 


IS 


21 


TC 


40 


A1 










41 


TC 


50 


A2 










51 


TC 


£C 


A3 










£1 


TC 


TC 


A4 










71 


TC 


EC 


A 5 











86 



ooonnnooonnooooooooooooooooooooonooooooooooooooooooooonooooonnnooooooooo 



A .4 TKC C IMENSIONAL I- E fiT TRANSFER ELEMENT DATA 




KST FCLLCW THE MATE MACPC CARD. 

THE TWC CIMENSIONAL HEAT TRANSFER ELEMENT IS CALLED 
ELMT02 ANC THUS THE ELEMENT TYPE NUMBER IN CCLLMN 10 
CF EACH MATERIAL SET NUMBEF CARD MUST BE 2 WHEN THIS 
ELEMENT IS REQUESTED. THE SECOND CARD GIVES GENERAL 
INFORMATION ANC IS PREPAREC AS FOLLOWS: 



. GENERAL INFORMATION CATA - FORMAT ( 5F10 .0 , 41 5 ) 



COLUMN 


< 


DESCRIPTION 


VARIABLE 


1 


TC 


ic 


CONDUCTIVITY IN X DIR. 
(IGNOREO IF TEMPERATURE 
DEPENDENT) 


D ( 1 ) 


11 


TC 


2 C 


CONDUCTIVITY IN Y DIR. 

( IGNOREC IF TEMPERATURE 
DEPENDENT) 


C < 2 ) 


21 


TC 


*0 


SPECIFIC HEAT (IGNOREO 
IF TEMPERATURE DEPEN- 
DENT) 


D (3 ) 


21 


TC 




DENSITY ( IGNORED IF 
TEMPERATURE DEPENDENT) 


D ( 4 ) 


41 


TC 


5 C 


HEAT GENERATION PER 
UNIT VOLUME (IGNORED 
IF TEMFEFATUFE DEPEN- 
DENT) 


C ( 5 ) 


51 


TC 


c c 


NUMBER CF INTEGRATION 
PCINTS FER DIRECTION 
(DEFAULT 4) 


NGP 


56 


TC 


6 C 


GEOMETRY TYPE 

.EQ.2 FCE AXISYMMETRY 

•NE • 2 FOR PLANE GEOMETRY 


KAT 


61 


TC 


65 


TOTAL NUMBER OF LINES 


NLBC 



WITH SPECIFIED BGUN- 
CARY CCNCITICNS IN 
ELEMENTS WITH THE SAME 
MATERIAL SET NUM8ER 

66 TC 7C CODE TO INCICATE IF ANY INOL 
CF THE MATERIAL PROPER- 
TIES IS TEMPERATURE DE- 
CEPENCENT. 

.EQ.O IF ALL MATERIAL 
PROPERTIES ARE CONSTANT 
•NE.O IF ANY OF THE MATE- 
RIAL PROPERTIES IS TEM- 
PERATURE DEPENDENT 



THE FCLLCWING CATA ARE RECUIREO IF ANY OF THE MATERIAL 
PROPERTIES IS TEMPERATURE DEPENDENT. THEIR ORDER IS 
CRUCIAL ANC TFCSE CARCS WHICH CORRESPOND TC CONSTANT 
MATERIAL FFCFEPTIES ML ST BE OMITTED. 



87 



nooonorfcnonooonnooonoooooonoooooononoonooooriooonooooooooooonrfcorinoooonoono 



MATERIAL FFCFEPTY COCE -FGFMAT(4I5> 

THIS C ARC IS ALLWAYS RECLIPED IF INOL.NE.O A N C MUST 
IMMEDIATELY FOLLOW THE GENERAL INFORMATION DATA CARD. 



COLUMNS 
1 TC 5 



6 TC 1C 
II TC 15 

16 TC 2C 



DESCRIPTION VARIABLE 

•EO.O CONSTANT CONDUCTI- IKX 
VI TY IN X DIR. 

•NE.O TEMP. CEP. KX 
•EC. 0 CONSTANT CONDUCTI- IKY 
VITY IN Y DIR. 

.NE.O TEMP. CEP. KY 
.EO.O CCSTANT HEAT CA IRCC 
PACITY (SPECIFIC HEAT * 

DENSITY ) 

.NE.O TEMP. CEP. HEAT 
CAPACITY 

.EO.O CONSTANT HEAT GE- IC 
NERATIGN PER UNIT VOL. 

.NE.O TEMP. CEP. Q 



CONDUCTIVITY IN' X DIRECTION TABLE 
THE TEMPERATURE DEPENDENT CONDUCTIVITY IN THE X 
CIRECTIGN (KX) IS CALCULATED BY LINEAR INTERPOLATION 
EETWEEN TV* C CONSECUTIVE ENTRIES IN THIS TABLE. THE 
FIRST CARC GF THIS SET CF CATA INDICATES THE NUMBER 
CF ENTRIES IN THE TABLE. IT MUST BE FOLLOWED BY THE 
SAME NUMBER GF CARDS V*ITH A PAIR OF VALUES (TEMPE- 
RATURE ANC CORRESPONDENT PROPERTY) IN EACH ONE. THEIR 
CRCER IS CRUCIAL. THE FAIRS MUST BE ORDERED ACCORDING 
TO THE INCREASING VALLE CF TEMPERATURE. 

OMIT IF 1KX.EC.0 



FIRST CARC -FORMAT! I 5 ) 

COLUMN J DESCRIPTION 

1 TC £ NUMBER CF ENTRIES OF 

THE TABLE TC BE INPUT 



VAPIAELE 

N 



TABLE CATA - FORMAT ( 2 F 1C . C ) 

FCR EACH FAIR CF VALUES A CARC MUST BE INPUT. THE 
TOTAL IS THE NUMBER SPECIFIED IN THE PRECEDING CARD, 
LOWEST TEMPERATURE IN THE FIRST CARD, SECOND LOWEST 
IN THE SECOND CARD, ETC. 



COLUMNS 
1 TC 1C 
11 TC 20 



DESCRIPTION 

TEMPERATURE 

CONDUCTIVITY 



VARIABLE 
XX ( I ) 

YY ( I ) 



.... CONDUCTIVITY IN Y DIRECTION TABLE 

PREPARED IN THE SAME WAY AS THE CODUCT IVITY IN THE 
X DIRECTION TABLE. 

CM IT IF IKY. EC. 0 



... . HEAT CAPACITY TABLE 

FREPAREC IN THE SAME WAY AS THE CONDUCTIVITY IN THE 
X CIRECT1CN TABLE. 

OMIT IF IRCC. EC. 0 



.... HEAT GENERATION PER UNIT VOLUME TABLE 

FREPAREC IN THE SAME WAY AS THE CONDUCTIVITY IN THE 
> CIFECTICN TABLE. 

CM IT IF IQ.EQ.C 



.... LINE BCUNCARY CONDITION DATA - FORMAT! 3 1 5 , 2F 1 C . 0 ) 
CMIT IF NLEC.EC.O. 

A CARC MUST INPUT FOR EACH LINE BOUNDARY CONDITION. 
IF THE SAME LINE IS SUBJECTED TO MORE THAN ONE TYPE 
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CF BOUNDARY CONDITION, A CARD MUST BE USED FCP EACH 
CNE CF TFESE TYPES. 

1FE TCTAl M.NEER OF CAPOS MUST BE NLBC. 

THE PROPERTY VALUE IS CEFINED: 

FOR CONVECTION - CONSTANT FEAT TRANSFER COEFFICIENT 
(IGNORED IF KBCOND (I ) .EC .4) 

FOR FLLX - FLUX PER UNIT AREA 

FOP RADIATION - PRODUCT CF EMISSIVITY BY STEFAN- 

BOLT ZM AN CONSTANT 

THE AMEIENT TEMPERATURE IS IGNORED FOR FLUX BOUNDARY 
CONDITION 



COLUMNS 
1 TO E 
6 TO 1C 



11 TO 15 



16 TO 21 
26 TO 25 



DESCRIPTION 
ELEMENT NUMBER 
CODE TO INDICATE BOUN- 
DARY CCNCITICN TYPE 
.50. 1 FLLX 

.EQ.2 CONVECTION (CONS- 
TANT COEFFICIENT) 

.50.3 RADIATION 

.50.4 CONVECTION (TEMP. 

DEP. CCEFF . ) 

CODE TO INDICATE LINE 
.50. 1 S=+l LINE 
.50.-1 $=-l LINE 

.50. 2 T=+l LINE 
.50.-2 T=-l LINE 
PROPERTY VALLE 
AMBIENT TEMPERATURE 



VARIABLE 
K5L (I ) 
KBCOND ( I ) 



KLINE ( I ) 



PRCFB(I,1) 
PROPB (1,2) 



89 



nnnnnOnnnnnnnnnnnnnnnnnnnnonnonnnnnnonnnnonnnnnnnnnnnnnnnnonnnnnnnnnnnnn 



A. 5 THREE CIMENSIONAL HEAT TRANSFER ELEMENT CATA 



/? 

/ . / 

/ . +T / 

11 . + 9 

/ 19 + / 19 

/ + / 

/ . + / 

✓ . + / 

• 21++++++ I ++++++ 

7 +...14.. 

. + 

• + 

20 . + 17 

15 + 

• + 

. R+ 






13 

/ 



MUST FCLLCW THE MATS MACRC CARD. 

THE THREE CIMENSIONAL HEAT TRANSFER ELEMENT IS CALLED 
ELMTG3 ANC THUS THE ELEMENT TYPE NUMBER IN COLUMN 10 
CF EACH MATERIAL SET NLMBEF CARD MUST BE 3 WHEN THIS 
ELEMENT IS REQUESTED. THE SECOND CARO GIVES GENERAL 
INFORM AT ICN ANC IS PREFAREC AS FOLLOWS J 



GENERAL INFORMATION DATA - FORMAT (6F10 .0 1 4 15 ) 



COLUMN 


c 


DESCRIPTION 


VARIABLE 


1 


TC 


1C 


CONDUCTIVITY IN X DIR. 
(IGNOREC IF TEMPERATURE 
CEPENDENT > 


D ( 1 ) 


11 


TO 


2 C 


CONDUCTIVITY IN Y DIP. 
(IGNOREC IF TEMPERATURE 
DEPENDENT) 


D ( 2 ) 


21 


TC 


2C 


CONDUCTIVITY IN Z DIR. 
(IGNOREC IF TEMPERATURE 
DEPENDENT) 


D ( 3 ) 


31 


TC 


40 


SPECIFIC HEAT (IGNORED 
IF TEMPERATURE CEPEN- 
DENT) 


D ( 4 ) 


41 


TC 


30 


CENSITY (IGNORED IF 
TEMPERATURE CEPENDENT) 


D ( 5 ) 


51 


TC 


60 


HEAT GENERATION PER 
UNIT VOLUME (IGNORED 
IF TEMPEFATUFE DEPEN- 
DENT) 


D ( 6 ) 


61 


TC 


65 


NUMBER CF INTEGRATION 
POINTS PER DIRECTION 
(DEFAULT 4) 


NGP 


66 


TC 


7 C 


GEOMETRY TYPE 

•EQ. 2 FOR AXISYMMETRY 

•NE. 2 FOR PLANE GEOMETRY 


KAT 


71 


TC 


75 


TOTAL NUMBER OF SURFACES 
WITH SPECIFIED BOUN- 
DARY CCNCITICNS IN 
ELEMENTS VITH THE SAME 
MATERIAL SET NUMBER 


NSBC 


76 


TC 


EC 


CODE TC INCICATE I F ' ANY 
CF THE MATERIAL PROPER- 
TIES IS TEMPERATURE DE- 
DEPENDENT. 

•EQ. 0 IF ALL MATERIAL 
PROPERTIES ARE CONSTANT 


INGL 
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•NE . 0 IF ANY OF THE MATE- 
RIAL PROPERTIES IS TEM- 
PERATURE DEPENDENT 



THE FOLLCWING CATA ARE REQUIRED IF ANY OF THE MATERIAL 
PROPERTIES IS TEMPERATURE CEPENDENT. THEIR ORCER IS 
CRUCIAL AM C THCSE CARDS WHICH CORRESPOND TC CONSTANT 
MATERIAL FFCPERTIES MUST BE OMITTED. 

MATERIAL PROPERTY CODE -F0PMAT(4I5) 

THIS CARC IS ALLWAYS REQUIRED IF INOL.NE.O ANC MUST 
IMMEDIATELY FOLLOW THE GENERAL INFORMAT I CN CATA CARC. 



COLUMNS 
1 TC 5 


6 


TC 


1C 


11 


TC 


15 


16 


TC 


20 



21 TC 25 



DESCRIPTION 

•EQ. 0 CONSTANT CONDUCTI- 
VITY IN X CIR. 

•NE.O TEMP. CEP. KX 

•EQ. 0 CONSTANT CONDUCTI- 
VITY IN Y DIR. 

.NE.O TEMP. CEP. KY 

• EQ* 0 CONSTANT CONDUCTI- 
VITY IN Z CIR. 

.NE.O TEMP. CEP. KZ 

•EQ .0 CCSTANT HEAT CA 
PACITY (SPECIFIC HEAT * 
DENSITY) 

•NE.O TEMP. CEP. HEAT 
CAPACITY 

•EQ. 0 CONSTANT HEAT GE- 
NERATION PEP UNIT VOL. 

.NE.O TEMP. CEP. Q 



VARIABLE 
I KX 



IKY 



I KZ 



IROC 



IQ 



CONDUCTIVITY IN X DIRECTION TABLE 
THE TEMPERATURE DEPENDENT CONDUCTIVITY IN THE X 
DIRECTION (KX) IS C ALCLLATEC BY LINEAR INTERPOLATION 
BETWEEN TWO CONSECUTIVE ENTRIES IN THIS TABLE. THE 
FIRST CARC OF THIS SET OF CATA INDICATES THE NUMBER 
OF ENTRIES IN THE TABLE. IT MUST BE FOLLOWED eY THE 
SAME NLMEER CF CARDS WITH A PAIR OF VALUES (TEMPE- 
RATURE ANC CCFPESPCNCENT PFOPERTY) IN EACH ONE. THEIR 
CRCER IS CRUCIAL. THE PAIRS MUST BE OROERED ACCORDING 
TC THE INCREASING VALUE CF TEMPERATURE. 

OMIT IF IKX.EC.O 



FIRST CARC -FOFMAT (15 ) 



COLUMNS 
1 TC 5 



DESCRIPTION 
NUMBER CF ENTRIES OF 
THE TABLE TC BE INPUT 



VARIABLE 

N 



TABLE CATA - FORMAT ( 2 F10. 0 ) 

FOR EACH PAIR OF VALUES A CARD MUST BE INPUT. THE 
TOTAL IS THE NUMBER SPECIFIED IN THE PRECEDING CARD, 
LOWEST TEMPERATURE IN THE FIRST CARD, SECOND LOWEST 
IN THE SECCNC CARD, ETC. 



COLUMNS 
1 TC 10 
11 TO 20 



DESCRIPTION 

TEMPERATURE 

CONDUCTIVITY 



VARIABLE 
X X ( I ) 

YY( I) 



CCNDUCTI V II Y IN Y DIRECTION 
PREPAREC IN THE SAME WAY AS 
X CIPECTICN TABLE. 

CM IT IF IKY. EC. 0 



TABLE 

THE CODUCTIVITY IN THE 



.... CONDUCTIVITY IN Z DIRECTION TABLE 

PREPARED IN THE SAME WAY AS THE CODUCTIVITY IN THE 
X CIRECTICN TABLE. 
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CM IT IF IKZ.EC.O 



... . FEAT CAPACITY TABLE 

PREPAREC IN THE SAME WAY AS THE CONDUCTIVITY IN THE 
X CIRECTICN TAELE. 

CMIT IF IPCC.EC.O 



.... FEAT GENERATION PER UNIT VOLUME TABLE 

FREPAREC IN THE SAME WAY AS THE CONDUCTIVITY IN THE 
X DIRECTION TAELE. 

CMIT IF IC.EQ.C 



SURFACE ECLNCARY CONDITION DATA - FORM AT (315 ,2F10.0 ) 
CMIT IF NSEC .EC .0. 

A CARD ML ST INPUT FOR EACH SURFACE BOUNCARY CONDITION. 
IF TFE SAME SURFACE IS SUBJECTED TO MORE THAN ONE TYPE 
OF BCUNCAF Y CCNDITICN't A CARD MUST BE USEC FCF EACH 
CNE OF THESE TYPES. 

THE TCTAL NLMEER OF CARCS MUST BE NSBC. 

THE PRC PERT Y VALUE IS CEFINED: 

FOR CONVECTION - CONSTANT HEAT TRANSFER COEFFICIENT 
( IGNCPEC IF KBCCNC ( I ) .EC .A ) 

FOR FLLX - FLUX PER UNIT AREA 

FCR RAC IATICN - PRODUCT OF EMISSIVITY BY STEFAN- 

BOLTZMAN CONSTANT 

THE AMBIENT TEMPERATURE IS IGNORED FOR FLLX BOUNDARY 
CONDITION 



COLUMNS 
1 TO 5 
6 TC 1C 



11 TC 15 



DESCRIPTION 
ELEMENT NUMBER 
CODE TC INDICATE BOUN- 
DARY CONCITICN TYPE 
•EQ.l FLUX 

‘ ' CONVECTION (COMS- 

CCEFFIC IENT ) 
RACIATION 
CONVECTION (TEMP. 
CCEFF. ) 

TC INDICATE SURF. 



VARIABLE 
KEL ( I ) 
KBCCNC (I ) 



16 TC 25 
26 TC -5 



• EQ. 2 
TANT 

• EQ. 3 

• EQ. 4 
DEP. 

CODE 

.EQ. 1 R=+ 1 

•EQ.-l F=- 1 
•EQ. 2 S=+ 1 

.E0.-2 S=- 1 

.EQ. 3 T=+l 
•EQ.-3 T=- 1 

PROPERTY VALUE 
AMBIENT TEMPERATURE 



SURF 

SURF 

SURF 

SURF 

SURF 

SURF 



KSURF ( I ) 



PROPB(Irl) 
PRCFB ( I » 2 ) 
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1 LINES WITH SPECIFIED BOUNDARY CONDITIONS 
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**END OF MACRO EXECUTION** 
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